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Project Summary 


Under the provisions of contract number NAS8-40194, which was entitled 
“TREETOPS Structural Dynamics and Controls Simulation System Upgrade", Oakwood 
College contracted to produce an upgrade to the existing TREETOPS suite of analysis 
tools. This suite includes the main simulation program, TREETOPS, two interactive 
preprocessors, TREESET and TREEFLX, an interactive post processor, TREEPLOT, and 
an adjunct program, TREESEL. 

The capability of the TREETOPS simulation package was extended by providing 
the ability to handle geometrically nonlinear problems through the use of a specially 
designed “User Controller". Software for computing the modal integrals from any finite 
element code (e.g., NASTRAN, ANSYS, etc.) was developed in two forms — Fortran source 
code and a MATLAB script file. 

A “Software Design Document", which provides descriptions of the argument lists 
and internal variables for each subroutine in the TREETOPS suite, was established. 
Additionally, installation guides for both DOS and UNIX platforms were developed. Finally, 
updated User’s Manuals, as well as a Theory Manual, were generated under this contract. 
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1 .0 Introduction 


The TREETOPS simulation package is a multibody dynamics tool capable of 
formulating and integrating the equations of motion for a wide variety of aerospace and 
mechanical systems. Special attention to implementation and integration of controllers 
(passive and/or active) into the state equations was paid in the initial formulation of the 
code — which sets it apart form other mechanical dynamics codes. 

Overtime, certain deficiencies and limitations in the TREETOPS code have been 
noted in the user community. Among these deficiencies was a lack of software control 
documents, limitations to geometrically linear response regimes, and a lack of user friendly 
modal integral tools. NASA contract number NAS8-40194 was established to address 
these concerns. 

This report will detail the activities undertaken in pursuit of the goals established 
for contract NAS8-40194. The Software Design Document, not included here for the sake 
of brevity, will be discussed in Section 2 of this report. The updates and additions to the 
TREETOPS Manuals are presented in Section 3. Section 4 details the simulation 
enhancements — which include the addition of geometric stiffening to the code as well as 
the development of user-friendly modal integral software tools. Section 5 will conclude this 
report. 
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2.0 Software Design Documentation 


The TREETOPS family of tools is composed of the following major software 
modules or programs: TREETOPS, TREEFLX, TREESET, TREESEL, TREEPLOT. 
TREETOPS is a time history simulation tool for multibody systems with active control 
elements. TREETOPS considers the total structure as in interconnected set of individual 
rigid or flexible bodies capable of undergoing large translations and rotations relative to 
one another. In addition to large rigid body configuration changes, the components of the 
system may simultaneously experience small elastic deformations. TREETOPS uses the 
assumed modes method to model the flexibility of the bodies. In this technique, the elastic 
deformation of each body is approximated through a linear combination of the products of 
shape functions and generalized time coordinates. The shape functions are typically a 
subset of normal modes derived from a finite element model of the component. 

TREEFLX is a program designed to compute flexible body input data for the 
components of TREETOPS from both COSMIC and MSC NASTRAN finite element 
packages. This data consists of the component shape functions, generalized mass and 
stiffness matrices, and modal integrals. The output of TREEFLX is the ‘problem. FLN' file. 

TREESET is an interactive preprocessor program to assist the user in entering data 
for the various TREETOPS programs. TREESET acts as an editor for adding new data 
to a file, deleting old data, or modifying existing data. The output of TREESET is the 
‘problem.INT’ file. 

TREESEL is also an interactive program to aid the user in the selection of a subset 
component modes which model the flexible elements of the system. Fewer modes reduces 
simulation run time and facilitates control system design. TREESEL is based on the 
modified component cost analysis method of model reduction. 

The TREEPLOT program reads the output data file ‘problem.OUT’ written by 
TREETOPS and plots or prints those variables selected by the user. It is a menu based 
post-processor based on keywords. Other versions of TREETOPS have been modified 
to output data in a MATLAB compatible file for plotting and analysis of simulation results, 
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thereby minimizing the use of TREEPLOT. 

Oakwood College constructed a functional hierarchical organization of these code 
modules with call and caller graphs. A detailed flow chart of the entire code was also 
developed. For each of the files which make up the TREETOPS family of programs, 
Oakwood College generated documentation which consists: filename, purpose, calling 
routine, called routines, common blocks, input parameters, output parameters, unused 
parameters, input global variables, output global variables, local variables. 

For the TREETOPS program, 296 files were documented. The documented main 
routine is shown below. 

1.131 MAIN 

INTERFACE : 

Called by : 

None 


Calls subprograms 

OPNTOP INPUT INITZE DYNAMI LINRZE 

OUTPUT RESAVE 


DATA : 

Commons included : 

DBP.F, DBSC.F, BUFFER. F, DBC.F 

Input Parameter Name List : None 

Output Parameter Name List : None 

Unused Parameter Name List : None 


Input Global Variables List 
matun integer*4 
dtout real*4 
tend real*4 

time real*4 


MAT file unit number 
Data inerval 
Simulation stop time 
Simulation time (sec) 
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tlin real*4 Time used in the 

linearization option, 

=0 if L is chosen, =TEND + 1 for Z and N 
toend real*4 End of the time interval for 

data 

plotting, =TEND 

tostrt real*4 Start time for data plotting 

(is 

usually equal to 0) 

Output Global Variables List : None 


DESCRIPTION : 


Found_in_f ile : 

maintops. for 


Purpose : 

This is the main program of treetops and it calls seven other 

major 

subprograms to perform the task of TREETOP software. The 
subprograms 
are: 

— OPNTOP : Open necessary files to support TREETOPS input 

and 


output data. 

— INPUT : Drive the data input function for TREETOPS by 
calling 

the appropriate input routines. 

— INITZE : Initialize the simulation and perform initial 
computation of simulation data. 

— DYNAMI : Execute the simulation dynamics by performing 

integration and propagation of the system state 

vector . 

— LINRZE : Generate the linearized system dynamics. 

— OUTPUT : Store simulation output data in two data files. 

Both 

files contain identical information, but one file is 
formated and the other is unformated. 

— RESAVE : Outputs event information at the end of the 
simulation. 
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For the TREEFLX program, 1 19 files were documented. The documentation for the 
main routine of TREEFLX is as follows. 

3.118 TREFLX 

INTERFACE : 

Called by : 

None 

Calls subprograms : 

DEDAMP ERRWRT EXTDAT MEMCHK NODCHK OPNDAF 
OPNFACE 

RDAGMT RDNAST RDRTMD TMCALC WRTFLX 
DATA : 

Commons included : 

PARI . F PAR2 . F 

Input Parameter Name List : None 

Output Parameter Name List : None 

Unused Parameter Name List : None 

Input Global Variables List : 
augy integer*4 (pnbb) 
bbid integer*4 
bbnu integer*4 
flxbnm integer*4 (pnbb) 
flxid integer*4 (pnbb) 
nmodbj integer*4 (pnbb) 
nnodbj integer*4 (pnbb) 
nuflbd integer *4 
dammth character* 2 (pnbb) 
redmth character* 1 (pnbb) 
error logical*4 
doesex logical*4 

Output Global Variables List : 
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bbid 

integer*4 

bbnu 

integer*4 

nmde 

integer* 4 

nnde 

integer*4 

error 

logical*4 

DESCRIPTION : 

Found in 

file : 


mainpr . for 


Purpose : 

This is the main calling program for the TREETOPS 
pre-processor 

TREEFLX. Major inputs are the TREESET interactive setup file 
problem. INT and the NASTRAN output data for each body. The 
major output is a problem. FLN file for TREETOPS. 


For the program TREESET, 59 files were reviewed and documented. The 
description of the main routine for TREESET is as follows. 


2.42 MAIN 

INTERFACE : 

Called by : 

None 


Calls subprograms : 

ADD CONVRT DELETE FCHECK GETF HELPS 

MODIFY OPNSET PRINT RBUF SAVE STOP 

DATA : 

Commons included : 

VARS . FOR BUFFER . FOR 

Input Parameter Name List : None 
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Output Parameter Name List : None 

Unused Parameter Name List : None 

Input Global Variables List : 
revno integer* 4 

revdat character*8 
ifunct integer*4 

Output Global Variables List : 
revno integer*4 
revdat character*8 

DESCRIPTION : 

Found in file : 
main. for 

Purpose : 

This is the main routine for TREETOPS pre-processor TREESET. 

The output of TREESET is the problem. INT file which is read 

by TREETOPS. 

The TREESEL program consists of 81 files. Again, the documented main routine 
of TREESEL is shown below. 

4.78 TRESEL 

INTERFACE : 

Called by : 

None 


Calls subprograms : 


ACTUA 

MPRT1 

CCOST DECMAT MEZER 

MTEQS2 SGET 

MEZERS 

MLTS2 

DATA: 

Commons 
PAR. F 

included : 

SEL.F TCA.F 
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Input Parameter Name List 


None 


Output Parameter Name List : None 

Unused Parameter Name List : None 

Input Global Variables List : 

gmf double precision (nacma,nacma) 

gdf double precision (nacma,nacma) 

gkf double precision (nacma, nacma) 

ba double precision (nacma, numax) 

pa double precision (numax, nacma) 

ra double precision (numax, nacma) 

idev integer*4 

iout integer*4 

iu integer*4 

nac integer*4 

nxva integer*4 

nxa integer*4 

inpop integer*4 

nx integer*4 

nr integer*4 

nu integer*4 

nxv integer*4 

nbody integer*4 

ad double precision (nxmax,nxmax) 

bd double precision (nxmax, numax) 

cd double precision (nrmax, nxmax) 

dd double precision (nrmax, numax) 

Output Global Variables List : 
gm double precision (nxvma, nxvma) 
gd double precision (nxvma, nxvma) 
gk double precision (nxvma, nxvma) 
gmf double precision (nacma, nacma) 

gdf double precision (nacma, nacma) 

gkf double precision (nacma, nacma) 

ba double precision (nacma, numax) 
pa double precision (numax, nacma) 
ra double precision (numax, nacma) 
nac integer*4 

nxva integer* 4 

nxa integer*4 
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inpop integer*4 

nr integer*4 

nxv integer*4 

runs integer *4 (nbmax) 

nine integer *4 (nbmax) 

nrs integer *4 (nbmax) 

nte integer *4 (nbmax) 

nre integer *4 (nbmax) 

nts integer *4 (nbmax) 

nbody integer*4 

idbod integer *4 (nbmax) 

rst logical*4 

net logical*4 

ad double precision (nxmax, nxmax) 

bd double precision (nxmax, numax) 

cd double precision (nrmax, nxmax) 

dd double precision (nrmax, numax) 

dataun integer*4 

DESCRIPTION : 

Found in file : 
treesel . f 

Purpose : 

This is the main routine for the pre-processor TREESEL 

designed to aid the user in the modal selection / model 

reduction process. 

The TREEPLOT program consists of 48 files which were documented by Oakwood 
College. The main routine for TREEPLOT is displayed below. 

5.23 MAIN 

INTERFACE : 

Called by : 

None 

Calls subprograms : . 
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DEFALT PLOT 


PRINT RESET SAVECM SETUP 


DATA : 

Commons included : 

FRAMES . FOR 

Input Parameter Name List 

Output Parameter Name List 

Unused Parameter Name List 

Input Global Variables List 
buf character*50 (1000) 

count integer*4 


None 

None 

None 


Output Global Variables List : None 


DESCRIPTION : 

Found in file : 
mainp. f 


Purpose : 

This is the main routine of Treeplot software. It resets data, 
display a menu and call the proper routine according to the 

user 

choice from the menu. 
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3.0 TREETOPS Manuals 


Revised TREETOPS theoretical and user manuals were developed by Oakwood 
College under this effort. The User’s Manual was prepared in a format approved by the 
NASA COTR to allow the user to exercise the TREETOPS package. Additionally, the 
User’s Manual modified to correct current deficiencies and changes. The manual 
presents detailed information necessary to use all options available in the program. 
Various examples which cover the main features were included as well as trouble shooting 
information and software compilation and installation notes. Example problems were 
worked in the manual and supplied on disk to validate the operation of TREETOPS. 

The current TREETOPS Theory Manual was reviewed and revised. Parts of the 
TREETOPS and CONTOPS User’s Manuals were included in the revised Theory Manual 
to provide background on how TREETOPS works. The notation within the Theory Manual 
was changed to be consistent. A section on geometric non-linearities was added to the 
manual. 
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4.0 Simulation Enhancements 


Under this contract, a general procedure for incorporating geometric stiffening into 
the TREETOPS simulation package was developed. An algorithm, FORTRAN program, 
and MATLAB M file were also developed to compute the modal integrals and other data 
required for the ‘problem.FLX’ file. These accomplishments are discussed in Sections 4.1 
and 4.2 through documents previously distributed to NASA and reprinted here in entirety. 
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4.1 Incorporation of Geometric Stiffening into TREETOPS 


4.1.1 Introduction. 

The purpose of this task was to develop a procedure for incorporating geometric 
stiffening into the TREETOPS simulation package. Although TREETOPS had a 
rudimentary capability of handling geometrically nonlinear beams that were spin-stiffened 
in previous versions, no general procedure applicable to arbitrary bodies was made 
available. 1 The procedure developed herein is capable of addressing geometrically 
nonlinear problems that arise in a typical multibody system setting — a setting in which the 
inter-body forces, as well as external forces, must be considered in the development of the 
geometric stiffness terms associated with each body of the system. 

A TREETOPS “user controller" is used to actually augment the underlying equations 
of motion of each body with the appropriate geometric stiffening terms. The modeling 
techniques (i.e., additions to the standard TREETOPS model of a structure) required to 
facilitate the geometric stiffening will be illustrated in several example problems. A brief 
introduction to the theoretical basis of the procedure, followed by some general guidelines 
for implementing the technique into a standard TREETOPS model will precede the 
example problems. 

As an introductory note of caution, it must be mentioned that the inclusion of 
geometric stiffness into a TREETOPS dynamic model (or any other dynamic simulation 


^t does appear that in Version 10, the subroutine GENMOD attempts to mimic 
the procedure of Banerjee and Dickens [1]; however, there are two problems with this 
implementation. First, it is not at all clear how to use the GENMOD code — it was never 
documented. Second, the method of Banerjee and Dickens is not sufficiently general 
to solve the general problem presented by multibody systems with externally applied 
loads [2]. 
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package) is generally a task that should be undertaken only by individuals with a relatively 
high degree of experience in both TREETOPS modeling and the finite element method. 
Whereas the required additions to the TREETOPS model ( <...>.int and <...>.flx files) are 
straightforward, the preparation of the required generalized geometric stiffness matrices 
is somewhat laborious and typically requires the use of a finite element package (e.g., 
NASTRAN) capable of solving geometrically nonlinear problems. 

4.1.2 Brief Theoretical Background. 

This section presents a synopsis of the underlying mechanics associated with the 
inclusion of geometric stiffness into TREETOPS. A complete derivation of the associated 
nonlinear equations can be found in Reference [1]. For the sake of clarity, the inclusion 
of geometric stiffness in the “standard” structural dynamics problem, i.e., 


M5t + Kx = F (1) 

will be used to develop the salient features of the new TREETOPS procedures. 2 The 
stiffness matrix in the above equation can be broken down into two parts, 

K = K L + K a (2) 

where K L is the usual linear stiffness matrix and K a is the geometric stiffness matrix. 

The geometric stiffness matrix is a function of the internal (stress) loads acting in 
the body (i.e., K a = K 0 (a)). For linear structures, the internal loads are a linear function 
of the applied loads. In other words, o = off^, f^). If we make the assumption that the 
inertial loads’ contribution to K a is due to the rigid body motions of the body in question, 
then the inertial loads can be developed directly from Newton’s 2nd Law and Euler’s 
equation for the lumped mass case, e.g., 


2 Reference [1] contains the analogous development for the fully nonlinear 
multibody equations of motion; however, the simpler development presented herein 
should allow the user sufficient depth of knowledge to address a wide range of 
geometrically nonlinear problems. 
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ff = -m * a* and f f = a - to * I* • w 


(3) 


where ff is the inertial force, tf represents the inertial torque, m *and n* are the mass and 
mass moment of inertia matrix, respectively, a* is the translational acceleration, and finally 
a and to are the angular acceleration and rate vectors. The “k” that is tagged to the terms 
in Eq.(3) indicates an association with the k m nodal body (or in simple configurations, the 
k" 1 particle) of the body. 

Since the acceleration term a k is being computed as if the body in question is rigid, 
then we can write 

**=« 0 + a*p*+wi'(i)xp* (4) 

where p* is the position vector from the body’s origin (“o") to point k. This is illustrated in 
Figure 1. 



Figure 1. Singe Body Position Vectors. 


Equations (3) and (4) can be manipulated to give the following forms of the inertial 
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forces and moments acting on node k. 



1 0 
0 1 
0 0 


0 0 Pl 

0 -p z 0 

1 P , -P x 


-P y 0 -Px -Px 
Px -P X 0 -Pjr 

0 -p z -p z 0 


Py 

Px 

0 


Pz 0 

0 Px 

Px Py 


o 

X 

o 

y 

o 

z 


OL 


w y 

«5 

“x“y 

U) CO 
x z 

CO CO 


( 5 ) 


and 


( 6 ) 


Several things should be noted about the form of the inertial loads given in Eqs.(5) 
and (6). First of all, the above equations are restricted to the lumped mass case (i.e., we 
have relied upon the fact that the mass and mass moment of inertia for nodal body k is 
explicitly known). This restriction shall be overcome shortly. The most important thing to 
note about Eqs.(5) and (6) is that the 3x12 coefficient matrix is a function of only the 
body’s geometry and inertial properties. All the required kinematic information necessary 
to compute the inertial loads is isolated into the trailing 12x1 vector. 

It can be shown [1] that when the body under consideration has a distributed mass 
matrix representation, the inertial forces and moments acting on a particular node can be 
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written as 


n>T<t> QRB v k T m y rUl) vUM) ytifj) y^m 


= -[»£*«» y*W) y«uj) yUkk) vutf) 0 ‘*> (/.*)] 


u x co 2 

(O y O)J 


U x (O y 

(0,0), 


where <t> QRB is the matrix of “geometric rigid body modes" and the following definitions are 
employed: 


y*(«,r) a 


y n(*.') a y f(*.r) * y*(r,«) 
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The y£(*,r)are y£, («,r)terms are defined analogously. 

The form of Eqs.(7) and (8) appear somewhat unwieldy, but in fact the procedures 
followed to compute the inertial loads are straightforward. It should also be noted that 
Eqs.(7) and (8), which are sufficiently general to address the distributed mass case, 
collapse to the simpler forms of Eqs.(5) and (6) when a lumped mass representation is 
assumed. 

For brevity, Eqs.(7) and (8) can be stacked and written as 



where f ^ is a 6^1 vector, F* is a 6x12 coefficient matrix, and A, is the 12x1 vector of 
kinematic variables. The above equation can be generalized to compute the inertial load 
vectors for the entire body (i.e., all nodes, not just node k) by essentially “stacking up” all 
of the nodal inertial loads. This process leads to the following equation, 


Inert ^ * A f 

i »*1 12 12*1 


( 12 ) 


where is the nxl vector of inertial loads acting on the body, F t is a n x l2 coefficient 
matrix, and A, is the usual 12x1 vector of kinematic variables. At this point, it is useful to 
think of the F, matrix as simply the container for 12 nxl inertial load vectors — each 
associated with a particular kinematic quantity contained in the A, vector. 

We are now in a position to expand the geometric stiffness matrix into a form that 
will be appropriate for introduction into the dynamic simulation. Since the K a matrix is 
function of the loads acting on the body (inertial and external), it can be expanded into the 
following form 

= K.(C * W = *«(U) * *.(U (13) 


By substituting Eq.(12) into Eq.(13), the geometric stiffness matrix that accounts for the 
inertial loading (i.e., the “motion stiffness”) can be written as 

F, V)*,V) - E K 0 {F,{i))A, (14) 

/ /“I 
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where the F,(i) vectors are simply columns of the F, matrix Note that there are 12 
geometric stiffness matrices associated with the inertial loads, namely K a ( F , (/)),/ = 1,12. 

The geometric stiffening due to external loading can be expanded in a manner 
similar to the inertial loading, i.e., 

/ \ 


*.('«») = “a 




E^(^(0)-uo 


(15) 




/ i-i 


The external loads in the above equation have been expanded into ne “locator" vectors 
(F mt [i)) and associated magnitudes (f mt (i)). It should be noted that there are ne 
geometric stiffness matrices, each associated with one of the ne external loads. For the 
multibody case, all inter-body forces acting through joints or hinges are considered 
externally applied forces as far as the geometric stiffening is concerned. 

By expanding the K a matrix into the forms represented by Eqs.(14) and (15), the 
time-varying nature of K a is isolated to the A,(i) = A,(i,t) and f mt (i) = f Mt (i,t) terms. This 
allows the 12 K a {F,(i)) and ne K 0 (F Mt [i)) matrices to be computed in advance of any 
simulation effort. Any finite element code capable of producing geometric stiffness 
matrices should be able to provide the required data. 

As is typical in the consideration of flexibility of multibody systems, we will consider 
the flexibility of any particular body to be expressed as a linear combination of shape 
vectors (e.g., modeshapes). Hence, a need arises for the calculation of generalized or 
modal stiffness matrices. For TREETOPS, the linear modal stiffness matrix is input 
through the <...>.flx file. The modal geometric stiffness terms, which are defined below, 
are read into the user controller wherein they are multiplied by the appropriate modal 
coordinates. Elements of the modal geometric stiffness matrices are defined as: 

K oi jt * *X, F ,(/))*, (16) 


and 


K Q2_ q j a <t£K 0 (F„,(/'))4> y (17) 

where <J> ; is the j* h modeshape vector. The modal force associated with the geometric 
stiffening terms can now be written as 
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12 


y = 


EKi,.A,(/,o+ ek; 


L /-1 


02 'mt 


L t m\ 


(18) 


where n is the vector of modal coordinates. The calculation of the modal force terms 
given in Eq.(18) is performed in the user controller. This will be demonstrated in the 
example problems presented later in this report. 

4.1.3 General Implementation Discussion. 

As a prerequisite to actually incorporating the geometrically nonlinear capabilities 
into TREETOPS for a particular dynamic system, the analyst should go through a mental 
screening process to determine the necessity of the nonlinear analysis. The essence of 
this screening test is the question, “Is a geometrically nonlinear response likely from this 
system?” In general, this question is very difficult to answer. For relatively stiff structures, 
geometrically nonlinear responses are uncommon. However, there are several situations 
in which geometric stiffening may play a significant role in the system response. For 
rotating structures, if the spin rate approaches (or exceeds) the first natural frequency of 
the structure, then the geometric stiffness terms are necessary in the analysis. For 
structures constructed from beams and plates, compressive loading in the body can 
significantly modify the stiffness of the structure. When the user finds himself in the 
unfortunate position of not knowing whether or not to expect a geometrically nonlinear 
response, the inclusion of the nonlinear stiffening terms is the safest course to follow 
despite the not insigficant effort involved. 

The initial step in building the geometrically nonlinear TREETOPS model is to 
decide which of the bodies in the system require the geometric stiffening terms. Consider 
a satellite/antenna system during a slewing manuever. Since the core body (satellite) is 
likely to be much stiffer than the antenna, only the antenna body would be considered 
likely to produce a geometrically nonlinear response. Hence, the geometric stiffening 
terms utilized in Eq.(18) need only be computed for the antenna body, not the core body. 

The overall procedure for integrating the geometric stiffening terms into TREETOPS 
is presented in Figure 2. For each body requiring geometric stiffening terms, the first step 
in the analysis is to compute the F,(i) inertial load vectors. A Matlab m-file (fimakeS.m) 
has been created which will produce these inertial load vectors. The initial portion of this 
Matlab function is shown below. 
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Step 1 . 


Step 2. 

Compute the 12 inertial 


Compute the K g 

load vectors (F,(/)) 


matrices associated 

associated with unit 


with the F,{ i ) vectors. 

values of the 12 A,(i) 

f i 



fimake5(idof_ 

tab, mode, phir, mass) 

% function fi = fimake5(idof_tab,nnode, phir.mass) 

% this m file function constructs the fi array (ndof X 12) from the 
% idof_tab , mass matrix , 

% and also the phir matrix. 

The inputs to the function are an idofjab array, which is an ndof x 2 array that maps each 
degree of freedom of the FEM model. The idofjtab contains, for each DOF, the associated 
node number and coordinate direction (x = 1, y = 2, ... , 0 Z = 6). For example, if the jth 
DOF in the FEM model represents the y direction deflection of the 205 th node then 
idofJab( j, 1 ) = 205 and idofJab(j,2) = 2. The scalar nnode is simply the highest node 
number in the model. Typically this is also the number of nodes in the model. The phir 
input to the fimake5 function is the ndof x 6 array of geometric rigid body modes. Mass 
is the ndof x ndof mass matrix. 

The next step in the analysis is to compute the individual K 0 ( F,(i )) and K g (F mt (/)) 
matrices. Note that there will be ne K 0 {F mt (i)) matrices, and up to 12 of the K e (F ( (/)) 
matrices that need to be computed. A typical procedure for obtaining these geometric 
stiffness matrices is to apply (in the FEM model) the F mt {i) or F,(i) load vectors to the 
structure and then let the FEM compute the resulting K g matrix. In MSC NASTRAN, a 
buckling solution or a nonlinear static solution (SOL66) can be run which forces the code 
to assemble and store the K 0 matrix in the model database. If one chooses to run SOL66, 
then the following DMAP alter can be added to a SOL63 restart run to obtain the K g 
matrix Note that the DMAP given here is valid for MSC Version 65, and will probably need 
adjustments to be valid in other releases. (Note: MSC’s technical support actually 
provided the DMAP, and will surely have the updated version available upon request.) 


23 



Figure 



2 . 

Overview of the Geometric Stiffening Data Preparation. 
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$ Restart from SOL 66 run 
READ 9 
ALTER 287 t 
$ Insert DMAP Alters 

CONO LNSTIF, ACON * SKIP THE ALTERS FOR SUPERELEMENTS 
PARAM //C,N,LT/V,N, NOLOOP/ V , Y , LOOP I D=0/ 1 $ READ LOOPID CARD 

COND LNSTIF, NOLOOP % NO LOOPID, GO TO LINEAR STIFFNESS 

DBFETCH /UGV, ESTNL, , ./SOLID/LOOPID//DBSET3 S DATA FROM SOL 66 RUN 
PARAML ESTNL//C,N,PRES////V,N,NONLK % CHECK PRESENCE OF ESTNL 
COND LNSTIF, NONLK $ NO ESTNL, GO TO LINEAR STIFFNESS 

TA1, MPT,ECTS,EPT,8GPDTS,SI LS,ETT,CSTMS,DIT/ 

EST,DESTNL,GEI,GPECT,ESTL/V,N,LUSETS/S,N,NOESTL/ 

S,N,NP/2/S,N,NOGENL/SEID/S,Y,LGDISP=-1/ 

V,N,NLAYERS=5 $ GENERATE TABLE FOR LINEAR ELEMENT 
COND JMPKGG.NOKGGX $ 

EMG ESTL,CSTMS,MPT,DIT,GEOM2S, , ,/KELM,KDICT, , , ,/ 

S, N , NOKGXX/0/0/0// ////// /////K6ROT S STIFFNESS FOR LINEAR ELEMENTS 
EMA GPECT,KDICT,KELM,BGPDTS,SILS,CSTMS/KBJJZ,/ S 

EMG ESTL,CSTMS,MPT,DIT, ,UGV,ETT,EDT/KDLEL,KDLDI , , , ,/ 

S,Y,NOO=1/0/0//NP S DIFFERENTIAL STIFFNESS FOR LINEAR ELEMENTS 
EMA GPECT.KDLDI ,KDLEL,BGPDTS,SILS,CSTMS/KDLGG,/-1/S 

ADD KBJJZ,KDLGG/KLTOT/ S 

EMG ESTNL, CSTMS, MPT, DIT,GEOM2S, , ,/KELMNL,KDICTNL, , , ,/ 

1/0/0//////////////V.Y.K6ROT S STIFFNESS FOR NONLINEAR ELEMENTS 
EMA GPECT , KD I CTNL , KELMNL , BGPDTS, SI LS, CSTMS/KB J JZNL ,/S 

ADD KLTOT.KBJ JZNL/KB1/ S 

EMG ESTNL, CSTMS, MPT, DIT, ,UGV,ETT,EDT/KDELM,KDDICT, , , ,/ 

1/0/0//NP/ $ DIFFERENTIAL STIFFNESS FOR NONLINEAR ELEMENTS 

EMA GPECT, KDDICT, KDELM, BGPDTS, SI LS.CSTMS/KBDJ J,/- 1/S 

$ 

S PRINT GEOMETRIC STIFFNESS MATRIX 
MATPRT KBDJJ// S 

OUTPUT* KBDJJ//- 1/1 1/1 S 

t 

ADD KB1 ,KBD J J/KJJZ/ S 

JUMP JMPKGG S SKIP LINEAR STIFFNESS GENERATION 

LABEL LNSTIF $ LABEL FOR LINEAR STIFFNESS GENERATION 
ENDALTER S 


If NASTRAN is not employed, then the user should consult the User’s Manual of the 
particular FEM package in question to determine the procedure for extracting the individual 
K g matrices. Note that K a also is known as the stress-stiffening matrix. 


Before leaving the discussion of the generation of the K Q matrices, a brief comment 
on the construction of the F mt {i) vector is in order. As an example, let us assume that 
there is an externally applied load acting in the z direction of node number 10 in the FEM 
model of the body. In this case, F mt {i) would simply be an ndof * 1 vector of zeros, with 
the exception of a 1 in the row corresponding to the z direction of node 10. In general, the 
F mt { i ) vector simply identifies (with the unit value 1) the location of the i* external load. 
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It should be recognized that often it is not necessary to compute all 12 of the 
K©(^/(')) matrices due to the natural constraints of the system. For example, if the 
problem at hand is a planar problem (e.g., x-y plane) then a®, to x , u> y , a x , and a y are all 
zero, as well as any term that contains them in the A, vector. In the case of a planar 
problem, since only 4 of the A, elements can be non-zero, only the 4 corresponding 
K 0 (F t {i)) matrices need be computed. Recognition of the system constraints reduces 
(in the planar case) the computation burden by 66%. 

Step 4 in Figure 2 is carried out by simply pre- and post-multiplying the individual K 0 
matrices by the matrix of modal vectors. In other words, the modal geometric stiffness 
matrices are simply computed as indicated in Eqs.(16) and (17). The modal geometric 
stiffness matrices should be saved to disk since they will be required in the user controller 
employed by the TREETOPS model. 

The next step in the modeling process of geometrically nonlinear problems is to 
build the <...>.int file. Only a few additions to the normal <...>.int file have to be made to 
accommodate the geometric stiffening procedure. The first is that nmode+1 “special” 
nodes must be added to the body, where nmode is the number of flexible modes used to 
model the flexibility of the body. A “jet” actuator must be located at each special node. 
The modal deflections of the special nodes (contained in the <...>.flx file) are designed 
such that the first special node can observe (and effect) only the rigid body modes, the 
second special node can observe only the 1st flexible mode — and no others, and so on. 
This concept is illustrated in the example problems. A set of translational accelerometers, 
rate gyros, and rotational accelerometers are attached to the first special node to measure 
the rigid body accelerations and angular rates of the body. A Matlab function called 
modall.m has been produced which will aid in the development of the <...>.flx file for the 
geometrically nonlinear problem. Provisions have been made in this function to 
incorporate the special nodes into the model. 

Before presenting the example problems, a brief discussion of the “user controller” 
which actually computes the geometrically nonlinear terms for the body is in order. 
Essentially, the purpose of the user controller is to carry out the operations indicated in 
Eq.(18). In order to do this, the user controller must be implemented to accept the 
a 0 , a, and to sensors as inputs and drive the “jet” modal actuators with the geometric 
modal forces as outputs. Internally, the a 0 , a, and to can be combined into the required A, 
components inside the user controller. For example, the input to the user controller from 
a sensor may be o) z , but the A,( 9) = to x . These ideas are consolidated in the following 
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example problems. 



4.1.4 Example Problems 

The following examples of geometrically stiffened problems illustrate the techniques 
discussed in this document. All supporting models and software are included in the 
accompanying disk. 


4.1. 4.1 Beam Spin-Up. 

This problem concerns itself with computing the response of a planar beam being 
smoothly spun up about a pinned hub. This is illustrated in Figure 3. The beam hub is 
constrained to follow the profile given below. 


6. = 





t i T. 


(19) 


where the maximum spin rate is = 4 rad/sec and the ramp time is T s - 1 5 seconds. 
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Figure 4. TREETOPS Model of Spin-Up Beam. 

Figure 4 illustrates the TREETOPS model. Since the flexibility of the beam was 
modeled by 3 modes, there are a total of 4 special nodes (nodes 4-7). The modal gains 
associated with the special nodes is given in the following table: 
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All rotational modal deflections (gains) are zero for the special nodes. From the table it 
can be seen that special node 4 (SN4) is capable of only sensing rigid body motions since 
it has zero gains in each of the flexible modeshapes. SN5 can sense (and effect) motions 
associated with mode 1 because of the unit 6x gain in this modeshape. However, SN5 
cannot sense modes 2 and 3. The other special nodes have similar modal sensing and 
effecting capabilities. 

At this point, the reason for “designing” the modal gains associated with the special 
points in the way illustrated by the previous table can now be made apparent. Since the 
user controller produces the modal forces associated with the geometric stiffening terms, 
the model needs to have a node (i.e., a special node) in which only the first mode is 
effected by any load applied to this node. Similarly, the model should possess nodes in 
which only the second (or third) mode is effected (driven) by an applied load. The reasons 
that the first special node has no flexible modal gains is two-fold. First, the motion 
stiffening terms need to have access to the rigid body accelerations and angular rates, 
thus this special node is the ideal location to place the appropriate sensors. Secondly, 
a node is required which can be used to drive only the rigid body motions of the structure. 
This is necessary because a force acting on SN5 drives both mode 1 and the rigid body 
modes. Since the geometric stiffening terms should not effect the rigid body motions, the 
rigid body modes must be undriven. This is accomplished by simply applying the negative 
of the load applied to the SN5 node to the SN4 node. Similarly, the negative of the SN6 
and SN7 applied geometric stiffening loads is also applied to the SN4 node. 

For the spinning beam problem, if Node 2 is considered the origin for the calculation 
of the F, vectors, then the only modal geometric stiffness matrix that is required is K* 01 , 
which is associated with the A, (9 ) = term. The a 2 term is not considered is due to the 
fact that the angular acceleration does not induce axial loads in the beam — hence the 
geometric stiffness matrix associated with a r is zero. 

The user controller subroutine, graphically illustrated below, is contained in its 
entirety in the accompanying software. 



•*2 

3 

-E*i 

i - 1 
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As a final note, it is important that the user include all five of the modal integrals that 
TREETOPS accepts in the <...>.flx file. These terms play an especially significant role 
in the system response when geometric stiffening is being considered. 

The effect of including the geometric stiffness terms in the analysis of the beam 
spin-up problem is dramatically illustrated in the following plot. 


Beam Spin-Up Example 











00 









-0.00 









X 

\ 



/ 




-0-20 

t 

c: 


\ 



/ 





\ 


y 

/ 




g - 0.40 

<B 

o 

A ^A 


N 


y 













.12 -O.oO 
Q 

A O A 



\ 

\ 








\ 






-O.oO 

4 AA 



\ 








\ 






-7 .00 | 

OA 

70 

5.1 

i 

JO 10. 

Time 

00 IS. 

[sec] 

i 

00 20. 


Includes K^Q No KjQ 



Figure 5. Effect of Geometric Stiffening on Tip Response. 


The solid curve in Figure 5 represents the case in which the K Q terms are included 
in the TREETOPS run. This curve has been shown to be correct [1]. The dashed curve 
indicates that neglecting the geometric stiffening terms in this analysis leads to an 
unbounded response — clearly an erroneous result. In fact, this erroneous result is what 
is expected from all multibody dynamics codes which have no accommodation for 
geometric stiffening. 
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4.1.4.2 


Beam Inside Rotating Hub. 


In this example, a flexible beam (which is identical to the beam used in the previous 
example) is cantilevered from the rim of a rotating rigid hub. The goal of this problem is 
to determine the critical spin rate, Q^for which the beam becomes unstable and ultimately 
buckles. The centrifugal effects of the spinning hub cause the compressive loads in the 
beam; however, the coriolis terms in the equations of motion play a role in determining O cr 
The analytical expression for the critical spin rate results in Q er = 2.186 rad/sec. 



Figure 6. Beam Inside Rotating Hub. 
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Since this example is a two body problem, it is important to consider the inter-body 
forces as external loads on the beam for the purposes of computing the geometric 
stiffening terms. This is illustrated in the Figure 7, which is simply a diagram of the primary 
inter-body force, F x . 

For this example, two K a matrices must therefore be computed. The first matrix 
must be computed to account for the rigid body spin rate of the hub (Cl), and the second 
is due to the F x “contact" load between the two bodies. 



The user controller for this case must implement the following equation for the 
modal geometric stiffening forces - 


& = * F x •Kq , 2 ) • n (20) 

where the F x term is contained in the FCONI array which is accessed in the user controller 
by adding the include file “..\inc\dbh.f into the USCC.FORfile. it should also be noted that 
the contact forces in the FCONI array are expressed in the inertial coordinate frame and 
should therefore be converted to the Body 2 fixed frame. 
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Figure 8 is a plot of the tip response of the beam for two different steady-state Q’s. 
The geometric stiffening terms are included in both cases. As can been seen, the beam 
has an unbounded response at Cl - 2.19 rad/sec, however, it is stable at Q = 2.18 rad/sec. 
This is in agreement with the theoretical result of Q w = 2.186 rad/sec. 

An interesting point to note is that there is no elastic response of the beam until 
t = 2 sec. This is due to the fact that a small “tap" is given to the beam at t = 2 sec. 
Without this “tap”, no response occurs because there is no initial imperfection (or 
disturbance) to “push" the beam into the solution bifurcation. This behavior is typical 
when dynamic solution techniques are applied to buckling problems. 


Baam Inslda Spinning Hub 
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Figure 8. Tip Response of Beam. 
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4.1. 4.3 


Plate Buckling. 


In this final example, the critical buckling load for a simply supported plate will be 
determined. The inertial loads do not play a role in the geometric stiffening terms; 
however, the running load across the top and bottom of the plate is obviously the prime 
driver for the calculation of the geometric stiffening terms. 

Figure 9 illustrates the buckling problem to be examined in TREETOPS. The 
flexibility of the plate is modeled with 10 modes; therefore, 11 special nodes are required. 
These nodes are collocated at the center of the plate, and have the designations of SN20 
through SN30. Node 32 is located slightly above the center of the plate, and is used for 
monitoring the plate deflection. 
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Figure 9. Plate Buckling Problem. 







The table shown below contains the modal gains for the special nodes. The gains 
associated with all other directions are zero. As can be seen, SN20 is the “rigid body” 
node since it has no gain in any flexible mode. SN21 through SN30 are designed to effect 
modes 1 through 10, respectively. 


Translational Modal Gains for Special Nodes | 


Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Node No. 

1 

2 

3 

4 

5 
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a 

D 
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SN30 
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0 

0 

SB 

0 

m 

0 

0 

1 1 


The user controller for this problem is straightforward. As indicated in the figure 
below, the only input to the controller is the running load, p. This load is generated by a 
Function Generator, and has a magnitude which is set in the <...>.int file. (The <...>.int is 
on the accompanying disk.) As is the usual case, the modal forces associated with the 
geometric stiffening term are the outputs of the user controller. 
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The user controller is responsible for computing the modal geometric stiffness terms 
according to the equation given below. 

^ = (p-0-n (21) 


Note that in this buckling problem, only a single geometric stiffness matrix is 
required — namely, the K 02 associated with a unit running load value (i.e., p = 1). 

Figure 10 is a plot of the transverse deflection of N32 for two different values of p. 
The disturbance in both cases is a simple pulse (or “tap”) delivered to N32 at t = 1 sec. 
As can be seen in the plot, the p = 0.170 Ib/in case exhibits an unbounded response, 
whereas the p = 0.165 Ib/in case is bounded. In other words, this TREETOPS study 
indicates that the critical buckling load is somewhere between p = 0.165 and 0.170 Ib/in. 
The NASTRAN buckling solution produces a critical value of p = 0.167 Ib/in for the plate 
model — which is supported by the TREETOPS solution. 
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Figure 10. Response of Simply Supported Plate to Running Loads. 


4.1.5 Summary. 

This document has demonstrated the successful integration of geometric stiffening 
terms into the TREETOPS multibody dynamics package. The procedure has been 
developed in such a manner as to accommodate not only motion stiffening (due to inertial 
loads), but also geometric stiffening due to external loads acting on the bodies. This 
allows the procedure to be applied to multibody systems in which the inter-body forces are 
significant. 

A set of tools have been developed (fimake5.m and modall.m) which ease the 
computational burden placed upon the user when calculating the inertial load vectors and 
the <...>.flx files. However, obtaining the individual geometric stiffness matrices remains 
the user’s responsibility. Fortunately, many of the popular FEM packages (e.g., 
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NASTRAN, ANSYS, etc.) are capable of producing the required geometric stiffness 
matrices. 

The current methodology for computing the modal geometric stiffening terms is 
carried out inside the “user controller” within TREETOPS. The inputs to the controller are 
the various rigid body accelerations and angular rate values, as well as the magnitudes 
of the external loads. The outputs of the controller are the modal forces associated with 
the geometric stiffness terms. The modal forces are applied through a series of “jet” 
actuators to a set of special nodes which possess modal gains “designed" in such a way 
as to make it possible to selectively drive the individual modes of the body in question. 
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4.2 Mass Integrals for Multibody Simulations 


4.2.1 Introduction 


This report will present a theoretical discussion on the calculation of mass 
or modal integrals common to the equations of motion of multibody systems. An 
algorithm and accompanying FORTRAN program will also be described to compute 
these terms from lumped or consistent mass matrices. 

4.2.2 Derivation of Mass Integrals 


Shown in Figure 11 is an undeformed generic flexible body. Let the 
continuous body be discretized into a system of p nodal bodies, with mass m, and 
inertia /, , and n total degrees of freedom; R is the position vector from the 
inertial frame O to the body fixed frame B; and p , is the body fixed vector locating 
the undeformed position of the nodal body I relative to B. Body deformation is 
defined through degrees of freedom at the nodes. Each node may have up to six 
degrees of freedom. The nodal degrees of freedom are placed in a state vector X. 
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x= 


DOF 1 
DOF 2 
DOF 3 


DOF n 


*i 

y i 



where x 1 is the x translational degree of freedom for node 1 and 6 V is the z 
rotational degree of freedom for node p. A degree of freedom map, available from 
most finite element programs, will be used to characterize the model. This map 
denotes for each degree of freedom a direction label, one through seven, and an 
associated node number. The degree of freedom type or direction is as follows : 

x= 1 y=2 z=3 

8 =4 e =5 0 =6 

x y z 


A value of 7 is used for generalized degrees of freedom in reduced or synthesized 
models. 

Figure 12 is the deformed flexible body. The translational and rotational 
deformation vectors of nodal body I are 5, and 0 f , respectively. These 
deformations are assumed to be small and compatible with linear strain- 
displacement relations. 

The deformed nodal body position relative to the inertial frame 0 is 
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Figure 11: Undeformed Flexible Body 



Figure 12: Deformed Flexible Body 







r, = R* p,*5, 


The inertial translational velocity of nodal body I is 

t t = R* u)x(p,*5 f )*6 f 

where to is the angular velocity of the body fixed frame B with respect to 0 and t, 
is the time rate of change of the deformation relative to B. 

The kinetic energy of the system can be written as 

pnocfa* * . 

The elastic deformations will now be described using the assumed modes 
technique. In this technique, the elastic deformations of the flexible body are 
approximated through a linear combination of the products of spatial shape 
functions and generalized time coordinates. The shape functions are typically a 
subset of normal modes derived from a discretized or finite element model of the 
component. 

For this derivation, there are at most n shape functions for the n degree of 
freedom discretized structure. The translational and rotational deformations for 
nodal body I are written as 

nmodw 

5 / = E V; 

ft mod** 

e /= E 4», ; n , 
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where 4> /; and ijj /y are the translational and rotational 3x1 vectors for the j-th 
shape function evaluated at the location of nodal body I, and n ; is the associated 
generalized time coordinate. This is essentially a coordinate transformation from 
the discretized displacements to a set of modal coordinates. The model is reduced 
in size by truncating the number of modes or shape functions used in the coordinate 
transformation. 

Let <t> be the matrix of eigenvectors of mode shapes. Each column of <t> 
corresponds to a mode and is associated with an eigenvalue. Each row of <t> 
corresponds to a model degree of freedom. 
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The following mass integral terms are defined as 

p nodes 

l”oy = E m i *i; 

/■ 1 


p nodes 

r i;= E ('”,P,*V l /%/) 

/■I 
p nodes 

k]k = E m i [ *</ */* 1 ‘ *// */* ) 
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V E m #(p/ r 4» v i -<t> v Pi r ) 

/■I 

pnodM 

*” 2 ;* = 5T m i < ^ij x ^ik 

/• 1 

where 1 is a 3x3 identity matrix. 

The computation of these terms is straight forward given the modes, DOF 
map, nodal geometry, and mass / inertia distribution at each node. The mass / 
inertia distribution is readily available for lumped mass matrices. However, most 
finite element programs employ consistent mass matrices. Consistent mass 
matrices are not diagonal as are lumped mass matrices. The mass is usually 
distributed using the same shape functions used in the finite element stiffness 
matrix. Therefore, an algorithm is desired to compute the mass integral terms for 
a consistent mass matrix. This is best arrived at through momentum arguments. 

The system linear and angular momentum vector can be written as 

P = MX 


where M is the mass matrix from the finite element model and X is the time 
derivative of the nodal degrees of freedom or nodal velocities. From the mass 
integral definitions, the total linear momentum vector for the system is 

nmodm* 

r oy n y 

7-1 

The total angular momentum vector for the system is 
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Now, the system momentum vector will be expressed in terms of modal coordinates. 

P = MX=M®r\ 

pnodmt 

P A = E P AI*P, XP T, 

/-I 

The system momentum vector for the j’th mode is simply 

P ; =M<t> ; n y = 

where P T'J and P A,J are the system linear and angular momentum vectors for node 
I and mode j. The system momentum vector can be written as 

nmodM 

p = E p j 

7-1 

The total system linear momentum vector is therefore 

nmodM pnodmt nmodM 

p t~ E E p tij = E Hyty 

7- 1 /-I 7-1 

Now define the product of the mass matrix and the mode shape vector as 

MOj=Aj 

The system momentum vector from the j’th mode can then be expressed as 

p j =a i^i 

where 
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A,jf\j\s the momentum vector at the i’th node from the j’th mode. The total linear 
momentum vector for the system is 


nmodm* pnodm* 

^T = ]C 5T A TIJ^J 

J - 1 /-I 


The first mass integral for the j-th mode is therefore seen to be 


pnodm* 

r 0 y = A tij 


From similar reasoning, the total angular momentum vector of the system is written 
as 


nmodm* pnodm* 

Pa* E E 

/•I /■ 1 


nmodm* pnodm* 

^"52 X) (A*/;* P/*^77;)ty 

y-i /-i 

The total angular momentum may also be written in terms of the second mass 
integral r l; in modal coordinates as 
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nmodtt 

p,* E r./l, 

The second mass integral is therefore seen to be 

pnodta* 

l"iy = 53 (A»//*P/ X ^ny) 

/- 1 

The remaining mass integral can be expressed in terms of A TIJ as 


pnod** 

^2. k ~ 53 A Tl jX<t>i K 

1 /-i 


^iy = 53 m i[f > i ^r/y^ ”^r/yP/ ) 

/■i 

Pn0dM / T- — V 

4 /(t = 53 m i {^tij $/* ^ ■ a t ,j <b tk ) 
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4.2.3 Software User's Guide 


A FORTRAN program based on the equations presented in Section 2.0 has 
been written to compute the mass integrals for discretized flexible bodies and 
output the TREETOPS .FLX file. This program is composed of the following files 
: massint8. for, input. for, flx_out.for, lib_flx.for, info.inc, 
inpt.inc, outpt.inc. The main program is in the file massints . for . It first 
calls the subroutine input to read the data, computes the generalized or modal 
matrices and mass integrals, and then call the subroutine f lx_out to write the 
.FLX file. The subroutine input is in the file input, for and flx_out is in 
flx_out. for . lib_f lx. for contains a library of matrix / vector manipulation 

routines, info.inc, inpt.inc, and outpt.inc are include files used to 
pass the necessary data between the various routines, info.inc contains the 
parameter statements for the maximum number of modes, nodes, and degrees of 
freedom. 

Upon running the program ( type massint depending on the compiler used), 
the user will be queried for the names of three files. The first file name is an ascii 
file containing information which describes the model size, degree of freedom map, 
discretized model node numbers, number of modes, and assorted output 
information. The contents of this file will be described in more detail in the next 
paragraph. The second file is also an ascii file which houses the discretized or 
finite element model mass, stiffness, and damping matrices. The matrices are read 
in input . for as follows : 
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open(unit=l, file = filel, status = 'old') 

cc Read mass, stiffness, and damping matrices in physical 
coordinates 

do i=l,ndof 

read(l,*) (mass_mat(i, j) , j=l,ndof) 
enddo 

do i=l,ndof 

read(l,*) (stiff_mat (i, j) , j=l,ndof) 
enddo 

if (idamp.eq. 1) then 
do i=l,ndof 

read(l,*) ( damp_mat ( i , j ) , j=l,ndof) 
enddo 
endif 

do i = l,ndof 

read(l,*) (phi ( i , j ) , j=l,nmode_sys) 
enddo 

close (1) 


The third file name is that of the output .FLX file. The name should be the 
root of the .INT file and include the .FLX suffix 

The model information file input data format is as follows : 


cc Read model size and FEM DOF map 
read (2,*) nnode , ndof , nmode_sys 
do i=l,ndof 

read(2,*) inod(i), idof_type(i) 
enddo 

cc Read nodal geometry matrix 

read(2,*) (r_offset (i) , i=l, 3) 
do i=l, nnode 

read (2,*) (node_mat(i, j) , j=l,3) 
enddo 
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cc Body id number 
read (2,*) bid 

cc Modal integral option number (0,1,2) 
read (2,*) modopt (l,bid) 

cc Body mass 

read (2,*) mass 

cc Damping matrix option 

cc 1 = read damping matrix in physical coordinates 
cc 2 = create damping matrix from diagonals of generalized 
cc mass and stiffness matrices 

read (2,*) idamp 

if(idamp .eq. 2) read (2,*) zeta 

cc Number of modes to output 
read ( 2 , * ) nmode_out 

cc Modal selection vector flag 
cc 1 = output modes in sequnetial order 
cc 2 = use selection array' 
read (2,*) nselect 

if (nselect. eq. 1) then 
do i = l,nmode_out 
nmdind(i) = i 
enddo 
else 

if (nmode_out . gt . nmode_sys) then 
write (6, *) 'trying to write out more modes than you 
read in! ' 

wri te( 6, *) 'review input deck - this programs is 
stopping' 

stop 
end if 

cc Read nmdout mode indices to be output 
cc These indices must be in ascending order 
do i = l,nmode_out 
read (2,*) nmdind(i) 
enddo 
end if 

cc Read in number of nodes to be output from discretized 
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model 


read ( 2 , * ) node_out 

cc Read in desired node numbers for modal output - phi and 
psi 

do i = l,node_out 

read ( 2 , * ) inod_out ( i ) 
enddo 

cc Read in special node option for geometric stiffening 
cc 1 = no special nodes for geometric stiffening 
cc 2 = output gains for special nodes 
read (2,*) spec_node 

The first line in this file reads in the number of nodes in the discretized 
model, nnode; the number of degrees of freedom in the discretized model, 
ndof ; and the number of modes in the file computed from the discretized model, 
nmode_sys . The second set of data is the degree of freedom map. For each 
degree of freedom, the software reads in an associated node number, inod(i) , 
and a degree of freedom type, idof_type (i) . The next line of data is the offset 
vector, r_of f set ( i) . This vector allows the user to compute the mass integrals 
about a point different from the origin used to compute the coordinates for the 
discretized nodes. The offset vector is labeled R in Figure 1 1 and translates the 
origin of the nodal geometry fromO to B. Note that TREETOPS assumes that the 
mass integrals are computed about the origin of the body fixed frame used to 
describe the nodal geometry in the .INT file. The next set of data are the 
coordinates for the nodes of the discretized model, node_mat ( i , j ) . Each row 
corresponds to the x, y, z components of the position for node i. The body 
identification number, bid, is the next line of data. After this line follows the mode 
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option number (0,1,2), modopt ( l , bid) . A value of 0 will direct the software to 
compute all mass integrals; a value of 1 will compute only the zero and first order 
terms; and a value of 2 will output only the zero order modal integrals. The body 
mass, mass, is input next to scale the alpha modal integrals. Following the mass 
is the damping matrix option flag idamp . A value of 1 directs the code to read in 
a physical damping matrix after the stiffness matrix. A value of 2 will enable the 
code to compute a diagonal generalized damping matrix from the product of twice 
the input value of zeta and the square root of the ratio of the generalized stiffness 
matrix diagonal and the generalized mass matrix diagonal. The next line of data is 
the value of zeta, zeta, if the value of idamp is 2. Otherwise, the next line is the 
number of modes to be output in the .FLX file, nmode_out. The user may select 
which modes to output through the modal selection vector. A value of 1 for 
nselect will direct the code to output the first nmode_out modes in the data file. 
A value of 2 will direct the code to read in the selection array nmdind(i)as an 
array of ascending output mode numbers, one per line. The next line of data is the 
number of discretized model nodes to be output, node_out. After this is the array 
of desired output node numbers, inod_out, in ascending order and one per line. 
Finally, the last line of data in the special nodes flag, spec_node, for use in 
incorporating geometric stiffness matrices into the model. A value of 1 for this flag 
will cause the code to output no special nodes. A value of 2 will direct the software 
to write out nmode + 1 special nodes. The first special node is the rigid body node 
and will have zero modal gains. The remaining special nodes are given a modal 
gain of 1 for x translation for one mode each. The special nodes are appended to 
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the end of the normal nodes in the .INT file. 

The user of this code may readily modify input, for for specific data 
bases. A set of example data files is supplied with this software to test it upon 
installation to a particular computer system. The model information file name is 
test. inp. The discretized data file name is test. mat. The output .FLX file is 
test . fix and matches the data computed using raodall . m in MATLAB. 
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5.0 Conclusions 


This report has described the enhancements and additions to the 
TREETOPS suite of dynamic analysis tools by Oakwood College. The Software 
Design Document, inclusion of geometric stiffness into TREETOPS, as well as the 
development of user-friendly modal integral tools have been presented. 

The Software Design Document provides a description of the argument lists 
and internal variables used in each subroutine of the TREETOPS suite. 
Geometrically nonlinear response problems can be addressed by incorporating a 
specially designed “User Controller* into the suite. Both Fortran and MATLAB 
forms of user-friendly modal integral tools were delivered under the provisions of 
contract NAS8-40194. 

Additionally, a new set of User and Theory Manuals were prepared under 
this contract. These items, along with the others listed, should provide the 
inexperienced user more confidence in the use of the TREETOPS package. The 
Software Design Document is intended to allow the advanced user a more rapid 
and complete immersion into the code. This will no doubt ease the burden of future 
analysts who attempt to provide further enhancements to the TREETOPS suite of 
codes. 
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